######################Simulations for posterior predictive p-values#############
#################################testing functions##############################

#generate the data and calculate the posterior predictive p-values (procedure A1--A3)#
ppppp = function(params)
{
####if params$bootm>0 then use bootstrap sd estimator, otherwise use asymptotic formula####
  if(params$bootm > 0)
    stu <- "boot"
  else
    stu <- "asym"
  expit = function(x) 
  {
    1/(1 + exp(-x))
  }
  
  Z.pred.alpha = function(x, alpha)
  {
    lpscore = as.vector(x%*%alpha)
    pscore  = expit(lpscore)
    n       = length(pscore)
    rbinom(n, 1, pscore)
  }
################################generate the data##############################
  dta <- dta_generator(params)
  t.t <- dta$t.t
  xat <- dta$xat
  t.f <- dta$t.f
  xaf <- dta$xaf
  y <- dta$y
#############################The observed test statistics#######################
  test.stat.obs = cal.tau(dta,params,stu)
  integrated_tx <- data.frame(t = t.t,x = xat)
  colnames(integrated_tx)[-1] <- paste("alpha",1:(ncol(integrated_tx) - 1),sep = "")
  f <- as.formula(paste("t ~ 0 +",paste(colnames(integrated_tx)[-1],collapse = "+")))
  
  est.alphas <- glm(dta$t.t~dta$xat+0,family = binomial(link = "logit"))
######posterior sampling of \alpha(the parameter of ps model) using the true model#######
  alphas.t = MCMClogit(f, data = integrated_tx, mcmc = params$m)
  
  integrated_tx <- data.frame(t = t.f,x = xaf)
  f <- as.formula(paste("t ~ 0 +",paste(colnames(integrated_tx)[-1],collapse = "+")))
######posterior sampling of \alpha(the parameter of ps model) using the wrong model#######
  alphas.f = MCMClogit(f, data = integrated_tx, mcmc = params$m)
  # for (i in 1:params$m) {
  #   test.stat.pred[,i] <- cal.stat.pred(alphas[i,])
  # }
  test.stat.pred <- matrix(NA,nrow = length(test.stat.obs),ncol = params$m)
  for (i in 1:params$m) {
######posterior sampling of \alpha(the parameter of ps model) using the wrong model#######
    dta$t.t <- Z.pred.alpha(xat, alphas.t[i,])
    dta$t.f <- Z.pred.alpha(xaf, alphas.f[i,])
    test.stat.pred[,i] <- cal.tau(dta,params,stu)
  }
  ppp = apply(cbind(test.stat.obs,test.stat.pred),1,FUN = function(stat.pred){mean(abs(stat.pred[-1])>abs(stat.pred[1]),na.rm = T)})  
  if(params$direct.p == T)
    ppp
  else
    test.stat.pred
}
#######################the normal approximation based testing###################
cal.pvalue.normal<- function(params){
################################generate the data###############################
  dta <- dta_generator(params)
  if(params$bootm!=0){
    params$normalize <- T
    tau <- cal.tau(dta,params,s = "boot")
  }
  else{
    params$normalize <- T
    tau <- cal.tau(dta,params)
  }
  if(params$twosided == F)
    1 - pnorm(tau)
  else
    2*(1-pnorm(abs(tau)))
}


###############################THE MAIN FUNCTION################################
##############use params$Method to determine which function above to use########
cal.pvalue <- function(params){
  pvalue.use <- switch(params$Method,"rand"=cal.pvalue.rand,"normal"=cal.pvalue.normal,"boot"=cal.pvalue.boot,"bootrand"=cal.pvalue.bootrand,"ppp"=ppppp)
  tryCatch(pvalue.use(params),error = function(e){NA})
}